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We discuss the physics of embolic stroke using a minimal model of emboli moving through the 
cerebral arteries. Our model of the blood flow network consists of a bifurcating tree, into which we 
introduce particles (emboli) that halt flow on reaching a node of similar size. Flow is weighted away 
from blocked arteries, inducing an effective interaction between emboli. We justify the form of the 
flow weighting using a steady flow (Poiseuille) analysis and a more complicated nonlinear analysis. 
We discuss free flowing and heavily congested limits and examine the transition from free flow to 
congestion using numerics. The correlation time is found to increase significantly at a critical value, 
and a finite size scaling is carried out. An order parameter for non-equilibrium critical behavior is 
identified as the overlap of blockages' flow shadows. Our work shows embolic stroke to be a feature 
of the cerebral blood flow network on the verge of a phase transition. 

PACS numbers; 87.19.xq, 87.10.Mn, 87.10.Rt, 64.60.ah 



I. INTRODUCTION 

Stroke is a common cause of death and disability, and 
is expected to become more prevalent as average life ex- 
pectancy increases [l| . A common cause of stroke is from 
embolisation, where small pieces of thrombus and plaque 
detach from the insides of diseased arteries and move 
through the vasculature to become lodged in arteries sup- 
plying the brain. Depending on the duration and position 
of the obstruction, embolic blockages can prove fatal. 

Doppler ultrasound embolus detection in individuals 
at high risk of stroke reveals that large numbers of em- 
boli typically enter the cerebral circulation in advance of 
a major stroke occurring Embolus detection and au- 
topsy studies of patients who have undergone cardiovas- 
cular surgery show that patients can experience several 
thousand emboli during typical heart surgery [3] (which 
carries a 3-10% stroke risk Q). The interplay between 
these emboli is governed by flow dynamics, and thus the 
onset of stroke can be considered as a complex process of 
arterial obstruction involving multiple blockages. How- 
ever, little is currently understood about the impact of 
multiple emboli leading to the onset of stroke. 

To gain a better understanding of the triggers lead- 
ing to the onset of stroke from multiple emboli, we de- 
veloped a probabilistic Monte Carlo approach to predict 
the severity of arterial obstruction [5| . Our model is non- 
trivial because blockages divert the flow of new emboli 
that enter the arteries (i.e. there is an effective interac- 
tion between emboli) . We predicted a very rapid change 
from free-flow with little occlusion to severely obstructed 
arteries as the size of the emboli increased, which we sus- 
pected was associated with a phase transition. In this ar- 
ticle, we examine the origins of this behavior more closely. 

This article is organized as follows. In section HIl we 
introduce our model and algorithm, and carry out a basic 
fluid dynamical analysis of the vascular tree to justify 
our flow weighting scheme. To gain more insight into 
the model, we determine the limiting properties of the 



(b) f 




No flow 



Flow 




Flow shadow 



FIG. 1: A representation of a small section of the cerebral vas- 
culature, (a) Arteries supplying the brain typically bifurcate. 
The number of branches in the vasculature is prohibitively 
large for a full fluid-dynamical analysis. For this reason, we 
develop a minimal model (b) where nodes on a bifurcating tree 
represent arteries which can become blocked. When emboli 
block nodes, blood flow is weighted according to the number 
of nodes at the bottom of the tree receiving unimpeded flow. 
The flow carries new emboli away from existing blockages, 
introducing an effective interaction between emboli. 



model in section lllll To better understand the motion 
of emboli through the tree, we examine space-time plots 
of blockages in a single realization of the model (section 
W} . We examine the time correlator and correlation time 
in Sec. IIVI The order parameter is examined in Sec. [V] 
Finally, we summarize our results in section IVTl 



II. MODEL 

A. Vasculature 

We treat the cerebral vasculature as a bifurcating tree, 
which is consistent with high resolution images showing 
that > 97% of all branches in the cerebral arteries are 
bifurcations @. In this study, we do not consider the 
largest vessels such as the circle of Willis (CoW). Since 
most emboli travel into the middle carotid artery (MCA) 
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after passing through the CoW, we consider our model 
to be a simphfied version of the cerebral arteries supplied 
by the MCA. Bifurcations can be characterized by the bi- 
furcation exponent, which is defined through the relation 
xj = 2x'^ where Xg and Xd are the respective radii of the 
source and daughter vessels In this study, we use a 
single bifurcation exponent of 7 at all levels of the tree. 
When the bifurcation exponent is greater than 2, the to- 
tal cross sectional area of the daughter vessels is greater 
than that of the source vessel. In many organisms, the 
bifurcation exponent is expected to be approximately 3 
Q. A schematic of our model vasculature is shown in 
Fig. [TJ We vary Mg (the total number of layers in our 
tree) so that the vessel sizes range from Xmax = 1mm to a 
few /im (xmin), which is similar to the dimensions of the 
arteries in the brain. Node diameters are Xm = a;„iin2™^^ 
with m the level of the tree (m = corresponds to the 
smallest arteries and m — Mi — 1 to the largest). The 
nodes are connected by segments representing arteries. 



B. Fluid dynamics and flow weighting 

Emboli are carried by the blood, and as will be dis- 
cussed in sec. Ill C[ the trajectory of emboli at a bifurca- 
tion is related to the relative flow in the branches. Since 
obstruction by emboli causes redirection of blood-flow, 
it is relevant to determine the ratio of flows in the net- 
work downstream from a bifurcation in the presence of 
blockages using a basic fluid dynamical analysis. 

For initial determination of how the flow is modified by 
blockages downstream of a bifurcation, it is convenient to 
assume that flow in segments of the tree obeys Poiseuille's 
law, 

AP = (^) , . R, (1) 

where /i is the blood viscosity, a is the radius of the artery, 
I is the length of the segment of the artery that is being 
considered, q is the volume flow through the vessel, R 
is the resistance of the vessel to flow, and A_P is the 
pressure drop over a vessel segment. By following all 
possible routes from the root node where the pressure is 
Pa to the capillary mesh where pressure is Py, a set of 
2^"^^ equations can be determined, 

N-l 

Rn,i{n)qn,i{n) — Pa — Pv (2) 

n=0 

which are supplemented by 2^^^ — 1 simultaneous equa- 
tions of the form 

qn,i — qn-l,2i-l — 2i = (3) 

derived from conservation of flow to give 2^ — 1 equations 
for the 2^ — 1 flows in the tree. The nomenclature qn,i 
indicates the ith node in the nth level. 
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FIG. 2: (a) Flow weighting on small symmetric trees where 
the drop in pressure obeys Poiseuille's equation, and the bi- 
furcation exponent 7 = 3. Flow resistances were assumed to 
change with tree level according to Eqn. 2] (b) Flow weight- 
ing on small symmetric trees where pressure differences have 
a nonlinear dependence on the flow according to Eqn. [5] and 
7 = 3. The flow weightings tend towards the linear weighting 
scheme with increasing system size. 

It remains to assign realistic values to the resistances. 
The lengths of arteries are proportional to their radii 
Z = [9| and the radii typically follow Murray's law 
a1_i — 2a^ with the bifurcation exponent normally 
7 = 3 [8] (note that in this section only, n = is the 
root node, with n increasing on entering the tree). Thus 
Rn oc a^"^ oc 2'^"/^. Replacing the coefficients, 

^" - 

ao is the radius of the root node. Following standard the- 
oretical practice, the ratio of segment length to segment 
radius is set as / = 20a, i.e. k = 20 p^ . 

It is straightforward to invert the simultaneous equa- 
tions to establish the proportion of total flow into a 
daughter vessel, / = q^/iqA + qs), the proportion of 
flow in direction A at a bifurcation as X (the ratio of 
free nodes in direction A to total free nodes) is changed 
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(the flow in daughter vessels A and B is qa and qb re- 
spectively) [15|]. The result is shown in figure Ufa). As 
the tree size grows, the relative resistance of the end ves- 
sels increases as compared to the root node, and the flow 
solution approaches the line f = X. For 7 < 3 the linear 
approximation is even more closely approached. 

Directly downstream from a bifurcation, Poiseuille flow 
is not fully formed, leading to an increased pressure drop 
which depends nonlinearly on the flow according to the 
formula [Till. 



AP = R 



Infi 



(5) 



here p the density of the fluid. We solve the set of equa- 
tions resulting from replacing Poiseuille's law for the non- 
linear pressure drop of equation [5] (using the results from 
the linear flow analysis as an initial solution) showing the 
results in Fig. [UJb) . Again, the line f — X is approached 
as the size of the tree in increased. The approach is less 
rapid than that found with the Poiseuille analysis, but 
for very large trees, a linear relation between the ratio of 
free arterioles to the ratio of flows is expected. 

The results of the basic fluid dynamical analysis indi- 
cate that a linear flow weighting scheme is a good approx- 
imation to the true fluid dynamics. The linear weighting 
approximation is strictly valid in the limit that the flow 
resistance in the capillary nodes is much higher than in 
any of the other nodes. It is also known from studies of 
vascular physiology that most pressure is dropped across 
the smallest arterioles (see e.g. Ref. [HI) giving further 
evidence for the validity of a linear weighting approxima- 
tion. 



C. Embolus trajectory at bifurcation 

The response of an embolus to flow at a bifurcation 
is complicated. In the limit that the size of an embolus 
approaches that of the blood cells, the ratio of the average 
number of emboli that travel into each of the branches 
must be identical to the ratio of flows. However, when 
emboli are large compared with the sizes of vessels, the 
proportion of particles flowin g in to the daughter branches 
does not have a simple form . 

Bushi et al. have performed experiments to establish 
the trajectories of emboli at a single bifurcation, a.nd in 
Fig. [3]we reproduce some of their data from Ref. [ij. Fig. 
Oshows r — Ea/ {Ea + Eb), which is the probability that 
an embolus travels into branch A ys f = QA/iqA + Qb), 
the proportion of flow in direction A at the bifurcation 
(the flow in daughter vessels A and B is qa and qb respec- 
tively, and the number of emboli traveling into branches 
A and B is Ea and Eb)- Data shown as solid points 
are taken from in vitro measurements by Bushi et al. 
using pulsatile flow , and the dotted lines are our pa- 
rameterizations of the data. We note that Bushi et al. 
measured Eb/Ea and qB/qA, which we have converted 
to p and / using Eb / Ea — l/r — 1 and qB/qA = 1// — 1- 
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FIG. 3: r vs / at a bifurcation. Data shown as solid points 
are taken from in vitro measurements by Bushi et al. [l3l |. 
Our parameterization of the experimental data via equation 
[6] is shown as the dotted line. 



The diameter of the pipes forming the bifurcation in the 
experiment by Bushi et al. were quite large (6mm for 
parent and daughter B, and 4mm for daughter A) which 
explains why their 'emboli' also needed to be large to see 
the any nonlinear trajectory weighting p^ . 

We suggest that the relation between probability that 
an embolus travels down branch A and the proportion of 
flow in branch A has the sigmoid-like form. 



tanh [&tanh"\2/ - 1) -I- c] -I- 1 



(6) 



where the parameters b and c depend on the embolus 
size [1^. This satisfies known limits, for example if 
there is no flow (/ = 0), then no emboli can travel 
down the branch (r = 0) and likewise if all the flow 
is in the branch (/ = 1), then there is a probability 
r — 1 that emboli flow into the branch. Furthermore, 
the behavior of small emboli is known. When 6=1, 
the linear weighting is recovered, corresponding to very 
small emboli (i.e. clots of a few blood cells) which gen- 
erally follow the flow of blood. If branches are sym- 
metric, then / and r are symmetric through the map 
f ■— f — 1 — f and r := r' = 1 — r, which is why we 
are confident of the sigmoid-like form. Thus for sym- 
metric bifurcations, the curve must go through the point 
r = 0.5, / = 0.5 and c = 0. Fits of equation [S] are 
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also shown in Fig. [3l The parameters of our fits are 
6= 1.47±0.07, c= -0.137± 0.025 for the 3.2mm emboU 
and b = 1.28 ±0.13, c = 0.170 ±0.097 for 1.6mm emboh. 
The small values of c indicate that embolus trajectory is 
more likely controlled by flow ratio than the asymmetry 
of the vessels. The 1.6mm data could be compatible with 
the line r = /, but the 3.2mm data are not. This indi- 
cates that for emboli of diameters of the order of half the 
vessel diameter, the linear weighting scheme is applica- 
ble. 

We examine the effects of the more complex weight- 
ing scheme on the crossover or phase transition by con- 
sidering the extreme limits of the parameterization with 
6 = l,c = (which corresponds to a linear trajectory 
weighting when emboli are small) and b oo, c = 0, the 
latter corresponding a strong interaction limit where em- 
boli always travel in the direction of greatest flow. Note 
that currently only symmetric bifurcations arc consid- 
ered. 

Since embolization affects flow, the paths of subsequent 
emboli entering the arterial tree are diverted away from 
existing blockages. This generates an effective interaction 
between emboli, implying that the onset of stroke should 
be considered as an interacting many-particle problem. 



D. Algorithm 

We introduce 'emboli' into the tree with size do, pro- 
duced with rate (s^^) and dissolving at a constant 
rate a = Ad/ At (mm/s) (consistent with spherical em- 
boli, where the volume dissolve rate is proportional to 
the surface area). We take At = Is. An essential com- 
ponent of the model is an effective interaction between 
incoming emboli and pre-existing blockages: at a bifur- 
cation, emboli are more likely to move in a direction with 
more flow. 

We use a Monte Carlo simulation to obtain numerical 
solutions of the model and determine average behavior. 
Our algorithm is iterative with the following steps; 

1. On any time-step, an embolus may be created in 
the root node of the tree with probability At/r. 
In this study, a single initial embolus size, do, was 
assumed. 

2. All end nodes are examined to determine if there 
is a blockage upstream. This determines the dis- 
tribution of blockages which are to be used in the 
following steps. 

3. The emboli move according to the following rules. 

(a) If the embolus is larger than the current node, 
it does not move. 

(b) If all arterioles downstream are blocked, the 
embolus may not move since there is no flow. 

(c) If the embolus is smaller than the current 
node, and there is flow downstream, 



i. The proportion of the tree that is blocked 
downstream from the embolus to the left, 
HA, is determined from the distribution 
of blockages in step (5] The proportion 
of unblocked branches in direction A is 
Xa^I- riA- 

ii. Likewise, the unblocked proportion of 
branches in direction B is Xb — 1 ~ riB- 

iii. The flow proportion in direction A is de- 
termined to be / = Xa/{Xa + Xb)- 

iv. The embolus then moves in direction A 
with probability determined from equa- 
tion [51 Otherwise it moves in direction 
B. 

4. All emboli dissolve, leading to a linear reduction in 
radius during each time step. Completely dissolved 
emboli are removed from the simulation. 



III. STEADY STATE AND LIMITING 
BEHAVIOR 

We can estimate the total distribution of emboli in 
the system using a steady state approximation. The size 
of an embolus a time t after introduction to the tree is 
d = do — ta, where c?q is the initial size of the embolus 
and t is less than to — do /a. Thus, the average number 
density of emboli with size d in the vascular system is 
N{d) = 9{do — d)/aT, where 6 denotes a step function. 
In clinical situations, there is typically a range of initial 
sizes centered about a single size. Simulations indicate 
that behavior of the model with a range of initial embolus 
sizes is qualitatively similar to that with a single size. 

The total number of blockages at a single level is = 

J2m/-, ""° N{x)dx, with the exception of the tree level 



with similar size to the initial embolus (2 

2(m+l)/7a 



<do < 



'a;min)where A'^m = J^^^^^^^^^ N{x)dx. Choosing 
a do that has the same size as a node, 



N„ 



2(m+l)/7 _ 2™/'^ / 



/ar. 



(7) 



assuming the condition do > 2*^'"+^)/'''a:niin is satisfied. If 
do < 2'"/'''a;inin there are no emboli of comparable size 
to the node and = 0. If do is not the same size as 
at least one of the nodes in the tree, a prefactor [do — 
a;min2"/'^]/[xmin2('"+i)/'' - a;mi„2'"/^] is associated with 
Nm for the largest set of nodes with blockages. 

Obstructions at different levels of the tree have differ- 
ent consequences for flow in the end arterioles. For ex- 
ample, an embolus blocking the root (source) node of the 
arterial tree prevents blood flow to all end nodes, whereas 
a blockage in one of the smallest nodes only affects flow 
in that node. In general, the proportion of exit (m = 0) 
nodes not receiving flow because of an obstruction in level 
misp^ = Nm2"'/Ni,,t. Where Mast = 2^'-^ is the to- 
tal number of exit nodes in the tree. We also determine 
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the node level corresponding to blockages by the largest 
emboli (Afbig) via x 

^^^2^bie/l = do,SO 



Mhig = 7log2((io/a;min)- 



(8) 



We have assumed that the initial embolus size is the same 
as the size of one of the nodes. 

In general, the blockage proportion can be determined 
from the n-point correlation functions, gi-.-n, 

A/big -1 A/big -1 A/big -1 

P = X/ ~ X/ 9nmPnPm + ^ gimnPlPmPn 



n=0 n^m,=0 

+ ... + (_l)A4.g-ig„^^ 



l^m^n—Q 
Afbig-1 

n P" (9) 



n=0 



To compute all the correlation functions analytically 
would require detailed knowledge of the probability for 
each configuration (an analogue of the partition func- 
tion). Since the form for that function is not known, we 
examine limiting cases to gain insight. In the low den- 
sity limit, there is plenty of space in the tree, and it is 
expected that emboli avoid each other (as though there 
is a strong interaction). In the very dense limit, space 
in the tree is limited and blockages from different emboli 
have a high probability of removing flow from the same 
arterioles, which is indicative of weak effective interac- 
tion. 

First, we examine the dilute limit where emboli avoid 
one another, so that all emboli block different parts of the 
tree, i.e. the emboli cause mutually exclusive halt in the 
flow as perceived from the exit nodes. The proportion of 
the tree that becomes obstructed considering the steady 
state conflguration is 



A/big -1 A/big- 

p= J2 P"'^ 



m=0 



m=0 



Mast 



(10) 



when J2fn=l 'A^m2"/Mast < 1, and p ^ 1 otherwise. 
The sum ends at Mbig — 1 since the node level corre- 
sponding to Aibig is blocked for an infinitesimally small 
length of time before a recently introduced embolus dis- 
solves and moves to the next layer Mbig — 1- Substituting 
the expressions for Mbig, the sum of the geometric series 
is, 



P ■ 



Xmin[2i/^-l]((rfo/a;mi„)(''+i)-l) 



aT2 



M, 



L(2(7+l)/7 - 1) 



(11) 



This result has potential for clinical significance, as it 
shows that for 7 > 2 a volume conserving break-up of 
emboli (do and r r/2) leads to a net re- 

duction in blocked arterioles. As previously discussed, 
7 « 3 in the vasculature. Our analysis indicates that clot- 
breaking therapy could be an important tool for treating 
acute stroke, provided it can be delivered quickly. Given 
that the model is in the early stages of development, we 



add the caveat that this result is merely indicative and 
should not be used directly for clinical applications. 

For a state where the positioning of emboli in all layers 
is independent (non-interacting problem) the total prob- 
ability of end arterioles being blocked is reduced by the 
overlap between co-existing blockages in different layers, 

p = po + (1 -Po)Pi + (1 - bo + (1 -Pq)pi])P2 + ■■■■ (12) 

The first term in Eq. [T^] is the probability of an embolus 
in level stopping flow to end arterioles. Emboli block- 
ing the next layer have a probability pi of blocking the 
remaining (1 — po) arterioles still receiving flow after the 
level blockages are taken into account and so on. Eq. 
1121 may be rewritten 



Afbig-l A/big -1 

P = Y Pn- Y PnPn 



(13) 



n— Q 

Afbig-l Afbig-l 
+ J2 PlPmPn + ■ ■ ■ + (-1)'^'^'"''^^ n 

There is a special value for the parameters in Eq. 
[T^ where no flow reaches end arterioles {p — 1 ) . This 
occurs when the level with the largest emboli (of size 
c?o) becomes fully saturated. This is because more em- 
boli block the largest nodes as shown by Eq. [7] and 
the largest nodes supply the most arterioles. Examin- 
ing the p = 1 case leads to an equation that relates the 
parameters of the system when it is just fully blocked, 
7VMbig-i2*^'''''~^/Mast = 1, which leads to a prediction of 
the embolus size corresponding to all end nodes blocked. 



2(7+i)/7«^2^^^-ix;^,i„/(2i/^ - 1) 



1/(7 + 1) 



(14) 



It is interesting to note that since dj^^ oc t this result 
leads to the same conclusion about clot busting as in 
the dilute (strong-interacting) limit, i.e. that breaking 
up emboli reduces blockage. We note that time depen- 
dent fluctuations about the average density have not been 
considered in our analysis. Such fluctuations slightly in- 
crease the embolus size required to completely block the 
tree, but become irrelevant if the correct scaling to an 
inflnite lattice is made. The point at which the tree be- 
comes fully blocked is similar to a percolation threshold, 
since there are no direct paths from the root node to the 
exit nodes (arterioles). 

We note that our flndings for the interacting and non- 
interacting limits indicate quite different behavior. For 
the low density (strong interaction) limit, emboli avoid 
each other. For the very dense (weak interaction) limit, 
space in the tree is limited and blockages from different 
emboli have a high probability of removing flow from the 
same arterioles. As shown by our analytic results, the 
two limits exhibit quite different dependence on embolus 
prevalence (which depends on embolization rate, dissolve 
rate and size) . Eq. [TO] for the dilute limit corresponds 
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FIG. 4: (Color online) Finite-size scaling showing the change 
in the form of the time autocorrelation function G(t) with 
varying gamma: (a) 7 = 2, (b) 7 = 3, (c) 7 = 4. Simulations 
are carried out using the linear trajectory weighting. The 
tree had 20 levels, and a = 3 x 10~* and r = 0.1. The lines 
are contours on G{t). On increasing 7, the timescale can be 
seen to increase, with the shape of the contours becoming less 
rounded. 



to all g — 0. Eq. [13] corresponds to all 5 = 1. Thus, 
we expect either a sharp transition or crossover behavior 
between the dilute and dense limits. 



IV. TIME CORRELATOR 

A key question in this article is whether there is 
crossover or transition between the dilute and dense 



FIG. 5: The curves of integrated correlation time become 
more sharply peaked on increasing 7. r = 0.1, a = 3 x 
10"". Curves are shown for various 7 and linear trajectory 
weighting. The tree had 20 levels, and a = 3 x 10~" and 
r = 0.1. dmax is the maximum in the integrated correlation 
time curve. 



limits. In a non-equilibrium system at criticality, the 
timescale associated with fluctuations would diverge (in 
this case, the length of time that an additional embolus 
added to the system contributes to an overall increase in 
the level of obstruction). However, if there is a crossover, 
then while there could be increase in the timescale of cor- 
relations, no divergence would be seen. The appropriate 
measure of such fluctuations is the time autocorrelation 
function, which is computed from the expression, 



G{t) = 



{{p{t) - p){p{0) - p)) 



{{p{Q)-pr) 

and the time correlation function according to, 

-, No 

r{t)^—J2{{p,{t)-p){pM-p)), 



(15) 



(16) 



where Pi{t) = 1 if end node i is blocked at time t and 
Pi{t) = otherwise. 

In a critical system, T{t) should have the form. 



r{t) oc — e 



1 



-t/c 



(17) 



where y is a constant, and C is the correlation time. If 
G{t) is positive at all times (i.e. there is no anticorre- 
lation) it is possible to estimate the correlation time by 
summing (integrating) over the autocorrelation function, 
C w At X]n=o G{nAt) (this measure is commonly known 
as the integrated correlation time). 

To examine the existence of a phase transition, it is 
necessary of investigate the behavior of the system as the 
size is increased (note that the finite size scaling of the 
order parameter is discussed in the next section). The 
correct way of increasing the effective size of the system 
is unconventional. It is necessary to make the size dif- 
ference between vessels smaller, which means that the 
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FIG. 6: (Color online) Change in time autocorrelation func- 
tion G{t) with varying gamma: (a) 7 = 2, (b) 7 = 3, (c) 7 = 4. 
Simulations are carried out using the step function trajectory 
weighting. The tree had 20 levels, and a = 3 x 10~* and 
r = 0.1. The time correlations become significantly longer 
and the peak in the contours becomes rapidly sharper as 7 
increases. 



emboli interact with more levels as they dissolve [17| . 
This can be done by increasing 7, since the difference in 
vessel sizes between levels is 2^^^^ — 1. Fig. |4] shows the 
change in time autocorrelation function G{t) with vary- 
ing 7, when the linear trajectory weighting is used. So 
that the simulations with different 7 can be compared, 
we use the dimensionless parameter do/'^50%- The size 
^50% corresponds to half of end nodes receiving no flow. 
We note that o?5q% is approximately proportional to 7. A 
distinct change in the shape of the contours can be seen. 



with the time scale increasing most rapidly towards the 
center of the graph as 7 is increased. This can be exam- 
ined in the context of the integrated correlation time. 

Figure [5] shows Monte Carlo results for the integrated 
correlation time. Increase in 7 has an effect on the corre- 
lation time consistent with a phase transition, with the 
hump associated with the maximum of the curve becom- 
ing more defined, and moving closer to the form consis- 
tent with a divergence, and thus a phase transition, for 
large 7. Time and memory constraints limit the maxi- 
mum size of the tree, so we have been unable to examine 
systems with 7 > 4. 

Finally, the time auto-correlators are examined in the 
context of the step function weighting of embolus trajec- 
tory, with the results shown in Fig. [5) Again an increase 
in the timescale is seen on increasing 7, with qualitatively 
similar changes in the form of the graphs to those seen 
for the linear weighting scheme. In the case of the step 
function weighting, there are regions of anti-correlation, 
so the integrated correlation time can not be used to ex- 
amine the timescale. Since there is evidence of a phase 
transition at large 7, it remains to find a plausible can- 
didate for the order parameter. 



V. SPACE-TIME PLOTS AND ORDER 
PARAMETER 

In order to examine how the model changes on go- 
ing through the transition, we followed the motion of 
emboli through the model arterial tree for specific real- 
izations of the embolization. Fig. [7] displays spacetime 
plots showing which capillaries (end nodes) do not receive 
flow during a single run when the linear flow weighting 
scheme was used. Panel (a) shows the dilute limit where 
emboli are able to avoid each other. Each 'wedge-like' 
pattern represents the lifetime of an embolus, which ini- 
tially blocks a large number of nodes, then dissolves and 
moves to block a smaller node. There is no clear spatial 
ordering between the emboli. In contrast (b) shows the 
dense limit, where space limitations force emboli to cast 
'flow shadows' on the same arterioles: emboli in different 
flow levels cause the same end nodes to lose flow. This 
indicates that overlapping blockages could be important 
in the transition. 

It is also of interest to see what differences the step 
function weighting scheme makes to the distribution of 
emboli in the vasculature. Fig. [5] shows spacetime plots 
for 6 — > 00 (step function trajectory weighting) with 
model parameters that are otherwise the same as for Fig. 
[T] The effect of the step function weighting is to strongly 
redirect the emboli away from blocked vessels. Thus the 
emboli become more evenly spaced in the tree, leading to 
a state that appears to be more ordered in space. Again, 
when large numbers of emboli are present in the network, 
the tree becomes sufficiently populated that emboli can 
no longer avoid each other, so the flow shadows from 
blockages overlap. 
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FIG. 7: Space-time plot showing fiow in end arterioles in 
the model with b = 1 (linear weighting model). The tree 
consists of Me = 18 levels with 7 = 2, do = 0.2mm and 
a = 10~^mms~^. Fully blocked regions are shown in black, 
and free flow in white. Panel (a) shows the dilute limit where 
emboli have available space to avoid each other. In contrast 
(b) shows the dense limit, where emboli are forced to occupy 
the same trees since there is no remaining space. This in- 
dicates that overlapped blockages are relevant to the phase 
transition. 

The spacetime plots in Figs. [7] and [8] indicate that the 
main change in the state of the model on increasing em- 
bolus size is the overlap between blockages. Therefore, 
we suggest that the following measure of the overlap be- 
tween flow shadows, /3, may act as an order parameter, 

P^{J2n,n,), (18) 
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FIG. 8: As Fig. [7] but with b ^ 00 (step function trajectory 
weighting). The emboli become more evenly spaced in the 
tree, leading to a state that appears to be more ordered in 
space. 



where = 1 if an embolus in level i stops flow in the 
end arterioles at level 0. 

If P were the order parameter, then the timescale asso- 
ciated with correlations would peak at the embolus size 
where the order parameter becomes finite. To investigate 
this, Fig. [9]shows Monte Carlo results revealing how the 
time autocorrelation function changes with initial embo- 
lus size (bottom panel) alongside the overlap parameter 
(top panel). The hump features in both P and the time 
autocorrelation function are related to prefactors intro- 
duced when the largest emboli do not have the same size 
as a node. The timescale associated with correlations is 
largest when do = 0.14mm, corresponding to the embolus 
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FIG. 9: (color online) Main panel: Time correlator as do 
is changed. Mi = 20, 7 = 1.736, t'^ = O.ls"^ and a = 
3 X 10~'*mm/s. Errors on individual points are approximately 
3% of the t = Q value (3 standard deviations). The timescale 
is significantly increased close to the transition. Top panel: 
The overlap between fiow shadows shows how the significant 
increase in timescale is associated with the onset of order. 
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FIG. 10: (a) Overlap parameter, (3 showing a second-order 
transition. The hump features are related to partially oc- 
cupied levels, (b) Percentage of blocked states and limiting 
behavior. The vertical lines indicate the sizes of nodes at 
each level of the tree. Mi = 18 with 7 = 2, = 0.1s~^ and 
a = 3 X IQ-^mm/s. (Inset: Mi = 20, 7 = 2.151, t'^ = Cls'^ 
and a = 10~^mm/s.) Where no error bars are visible, they 
are smaller than the markers. For these parameters, overlap 
between blockages suddenly increases at do = 0.18mm in a 
second order transition. The percentage of blocked end nodes 
is always non-zero. 



size where overlap of flow shadows becomes significant. 
This is a strong indicator that (3 is the order parameter. 

We also address whether the proportion of blockages 
could be the order parameter, rather than the overlap 
of blockages. Fig. fTUf a) shows numerical results for the 
overlap of blockages, (3. For these parameters, overlap 
between co-existing blockages suddenly increases at do — 
0.18mm. In contrast, panel (b) shows that the percentage 
of blocked end nodes is always non-zero and is unlikely to 
be the order parameter of the non-equilibrium transition. 
The limiting behaviors are also displayed on panel (b) for 
completeness and agree with the Monte Carlo results. 
Generally, the steady state approximation becomes more 
appropriate as system size increases. 

We complete this section by examining the blockage 
p and overlap parameter /J as 7 is varied to carry out 
the finite size scaling, as shown in Fig. [TlJ We re- 
iterate that the finite size scaling associated with this 
problem is unusual. Increasing 7 means that the emboli 
(which dissolve linearly with time) interact with increas- 
ing numbers of nodes. As 7 is increased, the change 
in the proportion of blocked end nodes on increasing 
do/d^oy^ is more abrupt. As we showed in the last sec- 
tion, the timescale associated with fluctuations becomes 
more sharply peaked when 7 increases. We also show the 
overlap parameter /? on a logarithmic scale. P changes 
over several orders of magnitude as do is changed. As 7 
is increased, the values of (3 decrease for c?o < ^50% and (3 
increases for initial embolus sizes do > ^^50%- Scaling of 
this type is consistent with a phase transition, and can 
be seen for both the step-function and linear trajectory 
weighting. 



VI. SUMMARY AND CONCLUSIONS 

We have discussed the statistical physics of cerebral 
embolization leading to stroke based on analysis of a min- 
imal model of particles moving through the cerebral ar- 
teries. We investigated whether the onset of stroke is 
associated with a gradual crossover or a phase transi- 
tion. Our model is non-trivial since the trajectories of 
emboli are weighted away from blocked arterioles, gener- 
ating an effective interaction. Examination of the time 
correlator showed that the timescale of fluctuations in- 
creases significantly at a critical value, consistent with a 
phase transition. Finite size scaling adds further weight 
to this view. We note that finite size scaling shows that 
the phase transition strictly occurs when the bifurcation 
exponent 7 = 00, since for reasons discussed earlier in 
this article the system has a finite size when 7 is finite. 
Thus crossover behavior is found when 7 is finite, but 
is controlled by the large 7 phase transition. The or- 
der parameter was identified as the overlap of blockage 
fiow shadows defined by /3 = QZi^jniUj) . Our model 
shows that the onset of stroke can be thought of as driven 
by non-equilibrium critical behavior. Further modeling 
incorporating increasingly realistic embolus parameters 
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FIG. 11: Blockage and overlap parameter as 7 is varied to carry out the finite size scaling. Simulations are carried out using 
both step function and linear trajectory weighting. The tree had 20 levels, a = 3 x 10~^ and r = 0.1. Some of the small do 
points are missing on the plots relating to the step function trajectory weighting, since /3 = 0. As 7 is increased, the values of 
/3 decrease for do < d5o% and (3 increases for initial embolus sizes do > d5o%. Scaling of this type is consistent with a phase 
transition. 



and anatomy promises to aid our understanding of em- 
bolic stroke, cerebrovascular disease and perfusion injury. 

We briefly summarize the differences and similarities 
between our work and the directed percolation problem 
[T^ . In the directed percolation problem, flow stops when 
the number of broken bonds reaches a critical value cor- 
responding to an absence of routes from source to exit. 
There are some crucial differences between the standard 
directed percolation problem, and the one described here. 
The first major difference is that the liquid is carrying 
particles which can break bonds and move the system 
closer to the percolation threshold. The second difference 
is that flow routes re-form after a characteristic time as 
emboli dissolve. The third crucial difference is that the 
properties of the tree downstream from the input influ- 
ence the direction in which bond-breaking particles are 
carried, so blockages are not randomly distributed. In 
this respect, our model represents an entirely new type 
of percolation problem. 

The existence of a phase transition in the flow could be 
of clinical relevance (although we add the caveat that our 
results should not be directly applied to medical treat- 
ment at the current time). Not only would an organized 



state have more blocked nodes than a randomly ordered 
positioning of emboli in the cerebral arteries, but also the 
increased time scales associated with fluctuations close 
to the transition would tend to increase the risk of brain 
damage during embolization. 

Further work will be carried out to improve the model. 
For example, all branchings are currently symmetric, 
whereas asymmetric branchings are more common in the 
vasculature, so an algorithm will be developed to con- 
struct more realistic trees. Also, we are currently carry- 
ing out detailed laboratory tests of the paths of emboli at 
single asymmetric bifurcations, and in a silica replica of 
the major cerebral arteries, and these results will guide 
realistic development of the computational model. 
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